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Bose-Einstein condensation (BEG) of a noninteracting Bose gas of N particles iir a two-dimensional 
box with Dirichlet boundary conditions is studied. Confirming previous work, we find that BEG 
occurs at finite N at low temperatures T without the occurrence of a phase transition. The 
conventionally-defined transition temperature Tb for an infinite 3D system is shown to correspond 
in a 2D system with finite to a crossover temperature between a slow and rapid iircrease iir the 
fractional boson occupation No/N of the grouird state with decreasiirg T. We further show that 
Tb ~ 1/ log N at fixed area per boson, so in the thermodynamic limit there is no significant BEG 
in 2D at finite T. Thus, paradoxically, BEG only occurs in 2D at finite N with no phase transition 
associated with it. Calculations of thermodynamic properties versus T and area A are presented, 
including Helmholtz free energy, entropy S, pressure p, ratio of p to the energy density U/A, heat 
capacity at constant volume (area) Gv and at constant pressure Cp, isothermal compressibility kt 
and thermal expansion coefficient Op, obtained using both the grand canonical ensemble (GCE) 
and canonical ensemble (CE) formalisms. The GCE formalism gives acceptable predictions for S, 
p, p/{U/A), Kt and Up at large N, T and A, but fails for smaller values of these three parame¬ 
ters for which BEG becomes significant, whereas the CE formalism gives accurate results for all 
thermodynamic properties of finite systems even at low T and/or A where BEG occurs. 

PACS numbers: 05.30.Ch, 05.30.Jp, 03.74.Hh 


I. INTRODUCTION 

Some of the thermodynamic properties of a noninter¬ 
acting three-dimensional (3D) Bose gas without internal 
degrees of freedom and confined in a cubic box are well 
known, where macroscopic Bose-Einstein condensation 
(BEG) occurs in the thermodynamic limit below a phase 
transition temperature The experimental work on 

BEG greatly increased after the initial discoveries in 1995 
of BEG in ultracold atomic gases confined to harmonic 
potential traps.These discoveries led to much ad¬ 
ditional theoretical work on BEG in ultracold gases in 
harmonic traps that are in general anisotropic along the 
three Cartesian axes, especially including the effects of 
boson interactions. 

Theoretical studies of BEG have been carried out for 
ID and 2D Bose gases,which are relevant to the above 
experiments on ultracold trapped atomic gases. Here we 
take the order parameter of a BEG phase transition to be 
the fraction Nq/N of the boson occupation of the ground 
state No- A BEG phase transition occurs if the depen¬ 
dences of Nq/N and associated thermodynamic proper¬ 
ties of the Bose gas on temperature T are nonanalytic 
at a temperature defined as the BEG transition temper¬ 
ature Te. In 1967, Hohenberg showed that BEG cannot 
occur in ID or 2D at finite temperature T in the thermo¬ 
dynamic limit for a homogeneous Bose gas.^® However, 
this result does not rule out BEG in the thermodynamic 
limit in inhomogenous Bose gases. Indeed, Bagnato and 
Kleppner showed in 1991 that a BEG phase transition 
occurs in ID and 2D noninteracting Bose gases in power- 
law-potential traps.Eurthermore, BEG can occur in a 
finite 2D system containing a finite number N of bosons 


at finite T. In such a system, BEG entails a smooth in¬ 
crease in the boson occupation of the ground state (and 
also excited states) with decreasing T with no BEG phase 
transition occurring. In particular, Ketterle and van 
Druten studied 3D and ID boson systems in a harmonic 
potential for finite In addition to confirming that a 
BEG phase transition occurs in the thermodynamic limit 
in ID, they found that a smooth increase of BEG occurs 
in ID systems with decreasing T and finite N, i.e., with¬ 
out a BEG phase transition occurring. 

Less well studied is BEG of noninteracting bosons con¬ 
fined to a 2D box with Dirichlet boundary conditions 
where the wave function of each boson is zero at the edges 
of the box. Ziff et al. presented results for the pressure 
versus volume p(U) and heat capacity at constant volume 
Cy{T) in the thermodynamic limit for dimensions 1 to 5.® 
Ingold and Lambrecht^® found that a BEG phase transi¬ 
tion does not occur in ID or 2D for a finite number N of 
bosons with N = 10® — 10^, even though BEG itself does 
occur at low T. Deng and Hui calculated Cv{T) for finite 
2D systems with 1 < < 10® and also found a smooth 

increase of BEG with decreasing T, with no evidence for a 
temperature-induced phase transition.®® With Dirichlet 
boundary conditions, one expects a nonzero pressure at 
T = 0 in a finite 2D Bose gas,®® whereas a zero pressure 
is obtained by setting the ground state energy to zero or 
by using the grand canonical ensemble (GCE) formalism 
instead of the canonical ensemble (CE) formalism. Such 
studies of bosons in a 2D box are not just of pedagogical 
interest, since over the past few years cold-atom traps 
have been constructed with 2D box-like potentials.®®“®^ 

Most theoretical studies of BEG in 2D have been on 
systems confined to harmonic traps. In these systems 
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pressure and volume are not relevant thermodynamic 
variables, and hence the associated thermodynamic prop¬ 
erties isothermal compressibility kt, coefficient of ther¬ 
mal expansion ap and heat capacity at constant pressure 
Cp are also not relevant. On the other hand, for bosons 
in a 2D box these thermodynamic variables and prop¬ 
erties are appropriate. Here we report a comprehensive 
study of the thermodynamics of the noninteracting Bose 
gas in a 2D box with finite N and Dirichlet boundary 
conditions using both the GCE and CE formalisms. We 
present studies of the various thermodynamic properties 
versus T, A and N, as well as of the populations of the 
ground and low-lying excited states. For finite N and A, 
all properties versus T and A must be analytic and finite 
as discussed above. Of special interest is how these prop¬ 
erties behave for T —> 0 and T —>■ oo, where in the former 
limit the properties should be physically acceptable and 
in latter limit should correspond to those of the (clas¬ 
sical) ideal gas. As is well known, the GCE formalism 
gives unphysically large fluctuations in at T < 

Hence in addition to extensive calculations of the ther¬ 
modynamics using the GCE formalism, we also present 
calculations performed using exact iterative expressions 
within the CE formalism. The latter calculations are 
reported for properties for which the GCE gives incor¬ 
rect results for the parameter regimes in which signifi¬ 
cant BEG occurs, and we compare the results of the two 
approaches. The CE formalism gives numerically exact 
and some analytically exact results for all properties. 

The calculation methods are discussed in Sec. H. In 
Sec. HI we calculate a quantity Te which is later shown 
to be the A^-dependent crossover temperature between 
weak and strong increases in the boson populations of 
the ground and low-lying excited states with decreasing 
T and/or A. We find that Te —?► 0 in the thermodynamic 
limit A^ —?► oo at fixed A/N. Hence in this limit signif¬ 
icant BEG condensation does not occur in the ground 
state or excited states at any finite T or A. Our re¬ 
sults obtained using the GCE formalism are presented 
in Sec. IV. The calculations of the fugacity and frac¬ 
tions of condensed bosons in the ground state and in a 
state in each of the first four energy levels are presented 
in Sec. IV A, of Cy in Sec. IV B, of the Helmholtz free 
energy F and entropy S in Sec. IV C, of the pressure p 
in Sec. IVD, and of the isothermal compressibility kt, 
thermal expansion coefficient ap and heat capacity at 
constant pressure Cp in Sec. IV E. 

Our calculations of properties for V = 1, 10, 100 
and 1000 within the CE formalism are described in 
Sec. V, which begin with calculations of the quantum 
state boson population statistics in Sec. V A. Calcula¬ 
tions of F and S are presented in Sec. VB, where we 
show that the finite values of S' at T —0 present in the 
GCE calculations for finite N are incorrect because exact 
CE calculations show that the entropy at T —0 is iden¬ 
tically zero for all finite N. The pressure is calculated for 
various N in Sec. V C, where we find nonzero values at 
T = 0 as expected, in contrast to the null values obtained 


using the GCE formalism, and quantify p(T = 0) versus 
N and A. The ratio of p to the energy density is found to 
be pf{UjA) = 1 exactly, in contrast to the strong deviati- 
ations that occur with the GCE formalism in the small T 
and/or A ranges in which BEC occurs. The calculations 
of kt, Qip and Cp are presented in Sec. VD, where we 
show that these properties are finite and positive at low 
T and/or A, in contrast to the divergent and/or negative 
values obtained from calculations of these properties us¬ 
ing the GCE formalism. A brief summary of our results 
is given in Sec. VI. 


II. METHODS 

A. Single-Boson Wave Functions and Energies 

The wavefunction ip oi a particle of mass m in a 2D 
square of side-length L and area A = with zero po¬ 
tential energy inside and infinite potential energy outside 
the square in the xy plane is given by the Schrbdinger 
equation as 


ip{x, y) = C sin(A:a;a::) sm{kyy), (la) 


where the edges of the square are at a; = 0, L and y = 
0, L, and C is the normalization constant. Continuity of 
the wave function requires the wave function to be zero on 
all edges of the square (Dirichlet boundary conditions on 
the wave function), yielding the quantized wave vectors 


k 


X 


rix'JT 


UyTT 


(lb) 


where rix-, riy = 1, 2, .... Thus the spatial distribution 
of the number density ~ ■0^ of the bosons inside the 
square is inhomogeneous. Periodic boundary conditions 
give incorrect energies for the low-energy quantum states 
which can modify their statistical and thermodynamic 
properties at low T and/or A for finite N. 

The kinetic energy A of a boson in the 2D box is quan¬ 
tized according to 


2m 


2mA 


(nl + nl). 


( 2 ) 


We put N noninteracting bosons into the square box and 
write A = N/{N/A). Thus we use the parameters N 
and A/N (the area per boson) as independent variables 
instead of N and A. Then Eq. (2) becomes 


E = 




( 3 ) 


The thermodynamic limit corresponds to V —^ oo at 
fixed A/N. We do not shift the energy scale so that 
the ground-state energy becomes zero, as usually done 
in calculations of the properties of the Bose gas, except 
when calculating the fugacity z using the GCE formal¬ 
ism where the energy shift does not affect the calculated 
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values of z (see Sec. IIB 3) or in some cases where we use 
the continuum representation for the high-energy quan¬ 
tum states. 

We introduce the dimensionless reduced parameter 7 
defined by 

^ A/N ^ 

^ (A/NU ^ ^ 


where {A/N)-^ is the value of A/N of the Bose gas at the 
Bose-Einstein crossover temperature Te to be defined in 
Sec. III. Thus 7 is the area per boson in units of the area 
per boson at T = Te, or equivalently the area normalized 
by the area Ae at T = Te- The reduced energy E is 
dehned using Eqs. (3) and (4) by 


E = 


E 

^bTe 



(5a) 


where fee is Boltzmann’s constant and the parameter a 
is defined as 


2mkBTE{A/N)E 


(5b) 


In Sec. Ill we show that a = a{N) and at large N one 
obtains a ^ In N. The reduced temperature t is defined 
as 



( 6 ) 


so from Eq. (5a) one obtains 


E _ E 
k^T t 



(7) 


Thus E/{kBT) is a function of the product yt, so we 
define the additional reduced parameter x as 


B. Grand Canonical Ensemble 

1. Distribution Function 


The Bose-Einstein distribution function for the average 
number of bosons with fugacity z in a quantum state with 
energy E at absolute temperature T is 

which in reduced parameters is 

/BE(E,t)= (11b) 

The fugacity is related to the chemical potential /r by 

z = ( 12 ) 


where the reduced chemical potential /r is defined as 


^bTe 


(13) 


Thus Eq. (11b) can be written 


fBE{E,t) 


1 

_ I ■ 


(14) 


An important consequence of Eq. (14) is that /be is 
invariant under a uniform shift of all energies, includ¬ 
ing both E and /i, by the same amount. Hence for 
all calculations involving the factor 
for convenience we set the ground state energy Eq with 
Ux = Uy = I to be zero, and then z is calculated with 
reference to this energy. Then Eq. (9) becomes 


E 


ajN) 

Nx 


{n1 + Uy 


2 ), 


(15) 


and the Bose-Einstein distribution function (11a) be¬ 
comes 


X = jt, 


and Eq. (7) becomes 


( 8 ) 


fBE{nx,ny,x,N) 


1 


2 ^exp[^(n2 -Png - 2)] - 1' 

(16) 


E 


a{N) 

Nx 


{nl+nl). 


(9) 


2. Density of Orbital States 


Using Eqs. (2) and (9), we also write 

E _ nl+nl 
knT g 

where 

Nx 2mkBT A 

^ a{N) 


( 10 a) 


( 10 b) 


Thus if we need to hold both T and A constant in a 
calculation such as in Sec. VB to obtain Eq. (70), one 
must hold g constant. 


In a continuum enumeration of the orbital quantum 
states, the number of these states A^states in a quadrant 
of a circle in n space is Ngtates = However, this 

includes states with (n^, = 0 , Uy ^ 0), [uy = 0 , nx ^ 0 ) 
and nx = ny = 0 for which the wave function in Eqs. (1) 
is zero. The number of such states is {2n + \)/2, where 
the states with nx = 0, ny > 0 and riy = 0 , nx > 0 are 
shared by two quadrants and the state with nx = 0 , = 

0 is shared by four quadrants. Correcting for these terms 
gives 

TT 1 

-^states — ^ ( 1 "^) 
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The density of orbital states in n space in the continuum 
representation is then 

= = (18) 
an 2 


3. Fugacity 

The fugacity z is determined from the requirement that 
the average number of bosons N in the system is equal to 
the sum of the average number of bosons in each quantum 
state, i.e., 


N= ^ fB^{nx,ny,x,N). (19) 


Using the energy expression in Eq. (15), this becomes 

„ 1 


E 


n^,ny = lZ lexp 


2 (Af) 


Nx 


{nl + nl- 2 ) 


- 1 


( 20 ) 


By specifying given values of N and of a{N) derived later 
in Sec. Ill, one can solve this equation for z{x, N). Since 
X = yt, one also has ^ = z{'-ft,N). On the other hand, 
when the fugacity appears by itself in an expression such 
as in Eq. (44) below where the fugacity is not multiply¬ 
ing the exponential of energy divided by /cbT, one must 
use the fugacity Zunshifted calculated from the unshifted 
energy levels by solving 


N = 


E 


-1 


na,,ny = l ^unshifted 


Nx VU 


- 1 


( 21 ) 


Comparison of Eqs. (20) and (21) gives 


■S^unshifted — ^ 6Xp 


2a{N) 

Nx 


( 22 ) 


so it is not necessary to do a separate calculation of 
-Zunshifted(a:, ^) if z{x,N) is already known. 

The boson occupation number Nq of the nondegenerate 
ground state with nx = ny = 1 is given by Eq. (16) as 


^0 = (23a) 

l — z 


The requirement that 0 < Nq < N gives the allowed 
range 


0 < z < 


N 

ivTi 


(23b) 


for X = oo and cc = 0, respectively. The fractional occu¬ 
pation of the ground state by the N bosons in the system 
is then 


No z 

N N{l-z)' 


(23c) 


In the continuum representation of the energy level dis¬ 
tribution, we use the density of states in Eq. (18) and 
Eq. (15) becomes 


E 



2 ), 


(24) 


where = ni+n^,. Thus the Bose-Einstein distribution 

X y 

function (16) becomes 


fBE{n,x,N) 


1 


z 1 exp 


a{n^ — 2 ) 


- 1 


(25) 


In 2D with finite N, Eq. (20) does not have an ana¬ 
lytic solution for z and therefore must be solved numeri¬ 
cally. Eurthermore, one cannot break up sums such as in 
Eq. (20) into the contribution of only the ground state 
plus an integral over the remainder such as is done for 
3D Bose gases in the thermodynamic limit because as 
we will see for the 2D Bose gas with finite N, in general 
significant BEC occurs in excited states in addition to 
the ground state. Therefore, one must include a signifi¬ 
cant number of states above the ground state in the sum 
and then carry out an integral over the remainder, and 
Eq. (20) for solving for z becomes 


V'‘ 


"max-^x 


E E 


rix — l riy — l 

pOO 


exp + nl-2)]-l 
V{n) 


z-iexp[^(n2-2)] -1 


dn. 


(26) 


After z(x, N) is determined, the fractional populations 
Ni/N, ..., N 4 /N of a quantum state in each of the first 
four excited energy levels, respectively, versus x and N 
are obtained using 

N- 1 

-j;^{x,N) =—fBEinxi,ny„x,N), (27) 

where fBE(nx,ny,x,N) is given in Eq. (16) and 
-I- n^. = 5, 8 , 10 and 13 for the first four excited en¬ 
ergy levels, respectively. To calculate these populations 
versus reduced temperature t at fixed reduced area per 
boson 7 or vice versa one uses the definition of x in Eq. ( 8 ) 
to replace x in the results by yt. 

The grand partition function Z is given by^ 

lnZ = -^ln(l-ze"-®’/'=^'^) , (28) 

i 

which for our system reads 


InZ = — 


E 


In 


1 — ze JVx ("x+"h 2 ) 


(29) 


where as discussed above we set the ground state en¬ 
ergy to zero in multiplicative factors (or 

^^-Ei/k-aT) appear in a calculation. The sum is 

evaluated by first calculating z{x,N) as described above 
and then replacing the sum to 00 by a sum from 1 to 
JT-max plus a numerical integral from rimax to 00 similar 
to the procedure used to arrive at Eq. (26). 
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C. Canonical Ensemble 

The partition function Q{N) within the canonical en¬ 
semble formalism for a system containing N noninteract¬ 
ing bosons is given by the recursion relation^®’^^ 

1 ^ 

= [Q(0) = 1], (30a) 

k=l 

where Qi{k) is the single-boson partition function for a 
modified temperature T/k given by 

Qi{k) = ^exp 


kE^ 


(30b) 


and the sum is over all quantum states i. Using Eq. (9) 
for Ei/k^T gives 


OO 

Qi(^) = exp 


^ fT’Tj — 1 


kd , 2 2 \ 


(30c) 



FIG. 1: (Color online) Parameter a in Eq. (5b) versus logjQ of 
the number N of bosons in the 2D system (open red circles). 
The empirical three-parameter fit of the data by Eqs. (35) is 
shown as the solid blue curve. 


This double sum can be expressed analytically as 


Qi{k) 



^ — ka/Nx 



(30d) 


where 0a{u,q) is a theta function that Mathematica de¬ 
notes as EllipticTheta. 

Another useful recursion relation within the canonical 
ensemble is for the average number hi (TV) of bosons oc¬ 
cupying a quantum state with energy Ei, given by^®’^^ 


ni{N) 


1 


N 

^exp 

fe=i 


kEj 

k^T 


Q{N-k). 


(31a) 


For the 2D Bose gas under consideration one obtains 


ni{N) = 


QiN) 


N 

E 


exp 


kcL / 9 2 \ 


Q{N-k). 


(31b) 

The computations in this paper were carried out using 
laptop or desktop computers and Mathematica software. 


III. CROSSOVER TEMPERATURE FOR 
BOSE-EINSTEIN CONDENSATION IN A 2D BOX 

The usual prescription for calculating statistical and 
thermodynamic properties of a three-dimensional (3D) 
Bose gas is to utilize integral representations of all sums 
over quantum states except possibly for the ground state. 
Thus the number of bosons in excited energy states above 
the nondegenerate ground state in n space is 

poo 

TVexc= / Vin)fBE{n,T,N)dn. (32) 

JO 


To obtain Te, one sets T = Te, z = 1 and TVexc = TV.®d8 
In 2D for n —>■ 0, the integrand becomes ^ so eval¬ 
uation of the integral at the lower limit n —> 0 gives 
^ lnn|„_>o = —OO. Thus the integral diverges logarith¬ 
mically for n —^ 0. However, this very slow divergence 
suggests that one should use a discrete sum over and 
Hy for the lowest energy levels instead of an integral over 
all n to determine Te. 

We confirmed that a finite Te can be obtained in 2D for 
finite TV if the integral over energies of the Bose-Einstein 
distribution function is replaced by a sum over the lowest 
energy levels with small Ux , ny and the integral formula¬ 
tion is used to sum over larger n as in Eq. (26). As 
discussed in the Introduction and will be demonstrated 
in Sec. IV A, this Te is a crossover temperature between 
weak and strong increases in TVq/TV with decreasing T in 
a 2D boson gas with finite TV, and is not a BEC transition 
temperature. 

We set z = I as in the 3D case and a; = I in Eq. (16) 
to obtain the Bose-Einstein distribution function for the 
calculation of Te 

/BE(n„%,T e) = ^ (33) 

for the sum and 

/BE(n,TE) = ^ (34) 

for the integral. Then we numerically solved Eq. (26) 
for the parameter a as a function of TV. The a val¬ 
ues reached constant values with increasing Umax by 
n-max ~ 500. A plot of a versus logj^Q TV for rimax = 500 
and log 20 TV = 0, 1, ..., 17 is shown in Fig. 1 and 
the values are given in Table I. One sees that a ap¬ 
proaches linearity in log^g TV at large TV. Therefore we 
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TABLE I: Parameter a in Eq. (5b) versus logj^Q of the num¬ 
ber N of bosons in the system obtained using Eq. (26) with 
Umax = 500 and z = x = 1. From calculations of a versus 
Umax we infer that the quoted values are accurate to « ±1 
in the last decimal place. The fitted values obtained from 
empirical Eqs. (35) are also shown, together with the percent 
deviations of the fit from the calculated a values. 


logio N 

a 

Q-fit 


0 

0.41539 

0.54089 

-30.21 

1 

1.11125 

0.96884 

12.82 

2 

2.17711 

2.08950 

4.02 

3 

3.50392 

3.47528 

0.82 

4 

4.98840 

4.99046 

-0.04 

5 

6.56312 

6.57845 

-0.23 

6 

8.19079 

8.21147 

-0.25 

7 

9.85178 

9.87427 

-0.23 

8 

11.5355 

11.5578 

-0.19 

9 

13.2356 

13.2563 

-0.16 

10 

14.9484 

14.9660 

- 0.12 

11 

16.6917 

16.6843 

0.04 

12 

18.4020 

18.4093 

-0.04 

13 

20.1444 

20.1397 

0.02 

14 

21.8735 

21.8745 

- 0.00 

15 

23.6312 

23.6128 

0.08 

16 

25.3885 

25.3541 

0.14 

17 

27.0824 

27.0980 

-0.06 


fitted the eighteen {logj^Q N, a{N)} data points by an em- 
prical three-parameter Fade approximant 




Q-fit 


^0 + P2(logioA^)' 

1 + Di logio ^ 


and obtained the fitting parameters 


FIG. 2: (Color online) (a) Logarithm to the base 10 of 1 — z 
versus x for a variety of boson number N values, where z is 
(35a) the fugacity, x = yt, 7 is the reduced area and t is the reduced 
temperature of the Bose gas. (b) Expanded plot of the data 
in (a) for a; = 0 to 2 . 


Po = 0.540886, P 2 = 0.948778, Di = 0.537570. 

(35b) 

In the limit of large N the fit gives gives aat = -3.2832-1- 
1.7649 logj^g iV. The fit is shown in Fig. 1 and the fit 
values and deviations of the fit from the data are shown 
in Table I. The magnitude of the deviation is seen to be 
^ 0.2% for 10"^ < iV < 10^^ with the deviation increasing 
to 30% for TV = 1. 

Equation (5b) gives Te for finite N to be 
^ 2mkB{A/N)Ea{N)' 

where the parameter a{N) diverges for TV —>■ 00 as shown 
above. Therefore from Eq. (36) and Eig. 1, Te decreases 
monotonically with increasing TV, and Te{N —>■ 00) = 0 
at fixed {A/N)e-, i.e., in the thermodynamic limit. 


IV. RESULTS: GRAND CANONICAL 
ENSEMBLE 

A. Fugacity and Fraction of Condensed Bosons 

The fugacities 2; calculated versus the parameter x = 
yt for finite systems with specific values of N obtained 
by solving Eq. (26) for 2 using Umax = 200-1000 and 
the respective a{N) values in Table I are presented as 
logio(l — z) versus x in Eig. 2(a), with expanded plots 
for X < 2 in Fig. 2(b). One sees for these finite systems 
that log]^g(l — z) shows noticeable increases with increas¬ 
ing X near x = 1 which become more pronounced as N 
increases. Since x = yt, if A = Ae (7 = 1), which is 
the value of A at Te, then the a; = 1 crossover occurs at 
T = Te {t= 1). From Eq. (23b), logio[linia:^o(l - z)] = 
— logio(A^o + 1) = ~logiQ(V + 1), which is verified in 
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0.8 0.9 1.0 1.1 1.2 

X 

FIG. 3: (Color online) (a) Fractional occupation Nq/N of the 
ground state versus x = for three N values, (b) Expanded 
plots of Nq/N versus x for N = 10^ to 10^® and x = 0.8 
to 1.2. 


Fig. 2. 

The fraction Nq /N of condensed bosons in the ground 
state versus x is now obtained from Eq. (23c) as shown in 
Fig. 3(a) for N = 10^, 10® and 10^®. With decreasing x, 
the data approach a linear behavior in x for x 1 as 
N increases, described by ^ = 1 — x. Expanded plots 
of data near x = 1 for = 10® to 10®® are shown in 
Fig. 3(b). The fraction of bosons in excited states for 
X < 1 is given by = 1 — = x (not shown). 

The approximately linear decrease of Nq/N versus x for 
X < 1 at large N is different from the behavior of the 
3D Bose gas in the thermodynamic limit which shows 
Nq/N = 1 - (r/TE)®/2. 

The data in Fig. 3 show that for these finite systems, 
BEG of bosons into the ground state occurs and increases 
smoothly and continuously with decreasing x. Hence 
there is no phase transition associated with BEG into 
the ground state (and low excited states, see below). Fur¬ 
thermore, according to Eq. (36) and the a{N) behavior 


in Fig. 1, Te —>■ 0 for —>■ oo and hence BEG does not 
occur in the thermodynamic limit. On the other hand, 
real systems do not contain an infinite number of bosons, 
and hence potentially observable BEG is expected to oc¬ 
cur for finite N, but with no BEG phase transition asso¬ 
ciated with it. 

The 3D Bose gas in the thermodynamic limit at T = 0 
has all N bosons in the ground state and none in the 
excited states. With increasing T, a macroscopic occu¬ 
pation of the ground state still occurs until the temper¬ 
ature (almost) reaches the BEG transition temperature, 
i.e., Nq = 0{N), but the occupation of any excited state 
i is Ni = 0(1). For T > Te, the occupations of all low- 
lying states states are about the same and of 0(1). In 
reduced dimensions with a finite number of bosons and 
no BEG phase transition, one expects this behavior to 
change. The ratio Ni/N of the number of bosons in an 
excited state with energy Ei is given by Eq. (27) as 


^{x,N) 


1 z 

^ + < - 2)] - z’ 


(37) 


where a{N) is given in Table I. For 1 = 0 with Ux = riy = 
1 one obtains Eq. (23c). Then from Eqs. (23c) and (37) 
one obtains the additional ratios 


Nq 


1 - Z 

{exp + n^. - 2)] - 1} -P (1 - z) ’ 


(38) 


where 1 — z is plotted versus x in Fig. 2. The Ni/NQ 
ratios versus x for -I- riy. = 5, 8, 10 and 13 and for 
N = 10®, 10® and 10®® are shown in Figs. 4(a), 4(b) 
and 4(c), respectively, along with the respective Nq/N 
versus x plots from Fig. 3. If the reduced volume per 
boson is fixed at 7 = 1, the parameter x is simply x = 
t = T/T'e- In that case, one sees that for finite N, the 
four excited states show Ni = 0{N) even when t <§C 1. 
Furthermore, with increasing N the occupation of the 
excited states occurs more rapidly with increasing t near 
t = 1 and the Ni/NQ values of the excited states at low 
temperatures decrease substantially. For the 3D Bose gas 
in the thermodynamic limit with t < 1 one has Ni = 0 {1) 
and hence Ni/NQ = 0 {\/N)^ which qualitatively differs 
from the 2D case especially for the smaller values of Al, 
whereas for t > 1 one has Ni/NQ k. 1 as in the 2D case 
at sufficiently high t. 


B. Internal Energy and Heat Capacity at Constant 
Volume 


The reduced internal energy is defined as 


U = 


U 

A:bT’e 


(39) 








1.0 





FIG. 4: (Color online) The ratio No/N versus x from Fig. 3 
and the ratios Ni /No of the number of bosons Ni occupying 
a state in each of the first four excited energy levels Ei {i = 
1 to 4) with respective quantum numbers ri^. + Uy. = 5, 8 , 
10 and 13 in Eq. (7) for (a) N = 10®, (b) N = 10® and (c) 
N = 10®®. For the ground state + riyg — 2. At fixed 
volume 7 = 1 one has x = t = T/Te and hence the plots are 
then versus t. 


The internal energy per boson divided by fen’T is then 
U U 


NksT Nt 


(40) 


1 E{nx,ny,N) 

= ^ 2 ^ - YY - fBE{z,x,nx,ny,N). 


Elx - 1 

Using Eqs. (9) and (16) this becomes 


E 




U _ 1 a{N) 

Nt~ N Nx ^^^^^^-iexp[^(n2+n2-2)] -1' 

(41) 

Similar to Eq. (26), we reformulate the sum as 


m 

1 a(iV) 
N Nx ' 


(42) 






max X 


E E 


nl + nl 


nx = l n, 


2-1 exp [-^{nl + n2 - 2)] - 1 


'D{n)r 


flmax z 1 exp 


a{n^—2) 

Nx 


-dn 


- 1 


where rimax was in the range 200 to 1000. For the large-a; 
region a: > 2, we replaced — 2 in the expression for the 
energy in the integral by and the density of states 'Diri) 
in n space given in Eq. (18) is replaced by 'D{n) = ^ 
so that the integral could be evaluated analytically in 
terms of polylogarithm functions. We utilized the same 
strategy for calculations of Z and other thermodynamic 
properties for x > 2 when integrals such as in Eq. (42) 
were to be evaluated. 

The heat capacity at constant volume (area) per boson 
is given in reduced units by 


Cy 

NkB 


{x,N) 


x§^{x,N) 


dx 


(43) 


where -^{x, N) is obtained from Eq. (42) and the partial 
derivative is obtained as the x derivative of a spline func¬ 
tion of a list of closely-spaced x^{x,N) versus x data. 
Shown in Fig. 5(a) are plots of Cy/NkB versus x = jt 
for N = 10® to 10®® and 0.001 < x < 20. For each N, 
the large-x (high temperature and/or large area) data 
approach unity as predicted by the classical equiparti- 
tion theorem for the two translational degrees of free¬ 
dom of a boson in 2D. Successively expanded plots of 
the low-a; data for x < 1.4 and a; < 0.1 are shown in 
Figs. 5(b) and 5(c), respectively. No sharp features are 
visible at a; = 1, which corresponds to t = T/Tb = 1 
and 7 = A/Ab = 1, as expected for BEC in these fi¬ 
nite systems where Tb is a crossover temperature rather 
than a phase transition temperature. From Fig. 5(c) one 
sees that Cy ^ a; at small x for N > 10®, whereas the 
data for iV = 10® show positive curvature. At sufficiently 
smaller x one expects an exponential dependence of Cy 
on X due to the energy gap between the ground and first 
excited energy levels. 
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FIG. 5: (Color online) (a) Normalized heat capacity per boson 
at constant volume (area) Cy/NkB versns x = "ft for several 
N values and 0.001 < x < 20. Expanded plots of the data in 
(a) are shown for (b) x < 1.4 and (c) x < 0.1. 


FIG. 6: (Color online) Normalized entropy per boson S/Nks 
versus x = "ft for (a) N = 10^ — 10^® and 0.001 < a; < 1 and 
(b) iV = 10® - 10“ and 0.001 < a: < 0.1. 


C. Helmholtz Free Energy and Entropy 


For each value of N shown in Fig. 5, Cy is approxi¬ 
mately proportional to x for a; ^ 1. It is of interest to 
know whether or not the entropy 5'(T = 0) = 0. To de¬ 
termine that one must first calculate the Helmholtz free 
energy F, given 

= ln2;unshifted(a:, A^) - j^\nZ{x,N), 

(44) 

where F = F/k-QT^, is the reduced free energy and the 
fugacity Zunshifted is for the actual unshifted energy levels 
as given in Eq. (22). Then using the definition F = 
U — TS, one obtains 


5 

NkB 


(x,iV) 


U F 

-—(x,N) — —(x,N), 
ivr ’ ^ NV ' 


(45) 


where -^{x, N) was calculated above from Eq. (42) as a 
prerequisite for obtaining Cy. 
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TABLE II: Reduced entropy and compression factor p — 
for a: ^ 0 versus log^Q of the boson number N, obtained 
from Eqs. (49) and (51), respectively. 


logio N 

ATfen (* 0) 

l^(^^O) 

0 

1.3863E-b00 

6.9315E-01 

1 

3.3510E-01 

2.3979E-01 

2 

5.6102E-02 

4.6151E-02 

3 

7.9083E-03 

6.9088E-03 

4 

1.0210E-03 

9.2104E-04 

5 

1.2513E-04 

1.1513E-04 

6 

1.4816E-05 

1.3816E-05 

7 

1.7118E-06 

1.6118E-06 

8 

1.9421E-07 

1.8421E-07 

9 

2.1723E-08 

2.0723E-08 

10 

2.4026E-09 

2.3026E-09 

11 

2.6328E-10 

2.5328E-10 

12 

2.8631E-11 

2.7631E-11 

13 

3.0933E-12 

2.9934E-12 

14 

3.3235E-13 

3.2236E-13 

15 

3.5649E-14 

3.4539E-14 

16 

3.6841E-15 

3.6841E-15 

17 

3.9144E-16 

3.9144E-16 


For X —>■ 0 only the ground state with = Uy = 1 is 
populated, and Eq. (29) then gives 

llnZ(x^0) = ^ln(N + l). (47) 

From Eqs. (23b) and (41) one obtains 




2 a 

Nx 


(48) 


This term turns out to cancel the identical term in 
Eq. (46). Using these results and Eq. (44), Eq. (45) gives 


S(x^Q) + 

NkB N ) 




(49) 


where the first term originated from In Zunshifted and the 
second came from In 2^(a: —>■ 0), i.e., both terms origi¬ 
nated from F. Shown in Table II is a list of values of 
versus N obtained using Eq. (49). The value for 
N = 10^ agrees with the value in Fig. 6(b). The value for 
N = 10® is just below our resolution limit in Fig. 6(b). 
The results demonstrate that within the GCE formalism, 
S(x —>■ 0) is nonzero for all finite N. However, this re¬ 
sult is not correct. We prove analytically using the CE 
formalism in Sec. VB that S(x —>■ 0) is identically zero 
for any finite value of N. 


Following calculation of F/Nt from Eq. (44), 
S(x, N)/NkB was obtained from Eq. (45) as shown for 
N = 10® to 10®® in Fig. 6 in the x ranges (a) 0.001 < 
x < \ and (b) 0.001 < x < 0.1. One sees that for 
N = 10® — 10®®, evidently S'(x —>■ 0) —^ 0, satisfying the 
third law of thermodynamics and also showing that the 
ground states for these N values are nondegerate. How¬ 
ever, in Fig. 6(b) one also sees that S{x —>■ 0) = const > 0 
for N = 10®, a surprising difference from the data 
for the larger N values. Therefore, the same result 
would presumably occur within the GCE formalism for 
N = 10® — 10®® at sufficiently small x with sufficiently 
high numerical resolution. One anticipates that there is 
only one way to put all N bosons into the nondegener¬ 
ate ground state with Ux = Uy = 1. Hence the entropy 
at T = 0 must be zero. The nonzero entropy calculated 
for N = 10® at X 0 therefore demonstrates that the 
GCE formalism can give incorrect predictions for thermo¬ 
dynamic properties for finite N in the quantum regime 
with small x, as found previously for the fluctuations in 
N at low T.25-30 y show analytically 

using the CE formalism that the entropy for a; —>■ 0 is 
identically zero for any finite IV. 

In order to determine the source of the nonzero entropy 
at a; = 0 within the GCE formalism, we examine the 
contributions from each term in Eq. (45) for the entropy 
at a; = 0. From Eqs. (22) and (23b), one has 

In ZunshiftedCa; ->■ 0) = + In ■ (46) 


D. Pressure 


Within the GCE, the pressure p is given by^ 

-^ = inZ. (50a) 

We define the reduced pressure p as 

Another reduced pressure is 

p(N, x) = ^ = ^ ln[2(iV, a;)]. (50c) 

This quantity for a gas is sometimes called the “compres¬ 
sion factor” in the literature (see, e.g.. Ref. 34). 

Shown in Fig. 7(a) are plots of p versus x = yt for 
several values of N. One sees that with increasing x, 
which corresponds to increasing area and/or temperature 
of the gas at fixed N, p approaches unity, as required 
since in these limits one must obtain the ideal gas law for 
which p = 1. An expanded plot of the data for x < 1 is 
shown in Fig. 7(b), where one sees that p{x —^ 0) = const 
for N = 10®. We can obtain an exact value for p(x —>■ 0) 
as follows. For T —>• 0, the ground state is populated by 
all N bosons. Equation (23b) gives the fugacity as z = 
N/(N+1). Then Eq. (29) gives \nZ{t 0) = ln(7V-bl). 
Using these results and Eq. (50c) one obtains 


p(N, X —>■ 0) 


ln(A^-b 1) 


N 


(51) 
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FIG. 7: (Color online) (a) Compression factor p = pV/NkBT 
versus the parameter x = 'yt for several N values and 0.001 < 
X < 60. (b) Expanded plot of the data in (a) for 0.001 < x < 
1. (c) The derivative dp/dx versus x for the same values of N 
as in (a) and (b). 


FIG. 8: (Color online) Ratio of the pressure p to the en¬ 
ergy density U/A versus x for boson numbers N = 10®, 10® 
and 10®. The exact value obtained from the CE formalism is 
unity for all values of x and for any finite N. 

A list of p{N, X —>■ 0) values versus N is given in Table II. 
For N = 1000, one obtains p{N,x —>■ 0) = 6.9088 x 
10“®, in agreement with Fig. 7(b). For the larger N 
values, p{N, x —>■ 0) is too small to resolve on the scale 
of the figure. The derivative dp/dx is plotted versus x in 
Fig. 7(c). For N = 10® to 10®®, the data show regions of 
X over which dp/dx = const and hence p is linear in x as 
also seen over the respective x ranges with less precision 
in Fig. 7(b). 

The reduced pressure p is calculated from the above 
values of p{x) at fixed N obtained from Eq. (50c) ac¬ 
cording to 

P{l,t) =-p{x). (52) 

7 

From dimensional considerations, one expects p oc U/A. 
As discussed later in Sec. V C, the exact analytic result 
for the noninteracting 2D Bose gas obtained from the 
CE formalism is p = U/A, or p/{U/A) = 1, for all N 
and X. Within the GCE formalism, this ratio is equal 
to p/{U/Nt), where p is given by Eq. (50c) and U/Nt 
by Eq. (40). The ratio p/{U/A) is plotted versus x in 
Fig. 8 for iV = 10®, 10® and 10®. One sees that the GCE 
formalism gives incorrect p/{U/A) ratios for N = 10® 
and 10®, with the deviation from unity increasing with 
decreasing N and x. Similar deviations must also occur 
for larger but finite N at lower x values than plotted. 
These deviations from unity again illustrate the failure of 
the GCE formalism to accurately predict thermodynamic 
properties for hnite N at small x where significant BEC 
occurs. 

Of particular interest for the thermodynamics are p 
versus 7 isotherms, p versus t isochores and 7 versus t 
isobars. These relationships are generated parametrically 
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from p{x) using Eq. (52) and the definition x = jt. For a 
p versus t isochore, one chooses a particular fixed value of 
the reduced area 7 and t is then obtained from x accord¬ 
ing to t{x) = x/ 7 . Similarly, for an isotherm, one chooses 
a particular value of t and 7 is obtained as 7 ( 3 ;) = x/t. 
In order to obtain a 7 versus t isobar, one chooses a par¬ 
ticular value of p. Then using 7 = x/t, Eq. (52) gives 


= 'Im 

Once t is determined for a given value of x, one uses 
7 ( 3 ;) = x/t(x) to obtain 7 for that value of t. 

Isochores of p versus t with 7 = 0.5, 1 and 1.5 are 
plotted in Figs. 9(a), 9(b) and 9(c) for N = 10^, 10® and 
10^®, respectively. One sees that below an TV-dependent 
temperature, the isochores for these three reduced areas 
for a given N are nearly the same. This means that in 
the respective t range, the pressure is nearly indepen¬ 
dent of area as will be seen explicitly in pressure versus 
area isotherms. At higher temperatures, the pressure de¬ 
creases with increasing reduced area. 

Isotherms of p versus 7 at fixed t = 0.5, 1 and 1.5 are 
shown in Fig. 10 for N = 10®, N = 10® and N = 10®. 
The plots for TV = 10® and TV = 10® show unphysical 
regions at low temperatures with positive slope, corre¬ 
sponding to a negative isothermal compressibility kt ac¬ 
cording to its definition for a 2D system given by 


1 

kt 


= -A 


dp{T,A,N) 

dA 


(54) 


Furthermore, the regions of 7 for which dp{'y)/d'y = 0 
for all three values of TV correspond to regions of infinite 
compressibility, which is unphysical for a finite noninter¬ 
acting Bose gas. 

Reduced area 7 versus reduced temperature t isobars 
withp = 0.1, 0.5 and 1 are shown in Figs. 11(a) and 11(b) 
for TV = 10® and 10®, respectively. The thermal expan¬ 
sion coefficient a-p is defined as 

In dimensionless reduced form this becomes 

dp = apTE = i . (55b) 

The isobars for TV = 10® in Fig. 11(a) exhibit unphysical 
regions of negative thermal expansion for p < 1 and small 
7 that are not apparent in the isobars for TV = 10® and 
larger TV. 

The above unphysical predictions of the GCE formal¬ 
ism for the thermodynamic properties at low values of TV, 
t and/or 7 at which significant BEC occurs are rectified 
in Sec. V below when we consider the predictions of the 
CE formalism for the same properties. 



t 


LU 



LU 

> 

d. 




FIG. 9: (Color online) Reduced pressure p = pyB/iksTs) 
versus reduced temperature t = T/T-e. for reduced area 
7 = v/ve = 0.5, 1 and 1.5 for boson numbers (a) N = 10®, 
(b) N = 10® and (c) N — 10®®. Note the different scales for 
the ordinates in (a), (b) and (c). 
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FIG. 10; (Color online) Reduced pressure p = pv^/k^T^ 
versus reduced area 7 = u/ue for reduced temperatures t = 
T /Te = 0.5, 1 and 2 and boson numbers (a) N = 10^, (b) N = 
10® and (c) iV = 10®. 



t 



t 

FIG. 11: (Golor online) Reduced area 7 versus reduced tem¬ 
perature t for reduced pressures pred = P = 0.1, 0.5 and 1 and 
boson numbers (a) N = 10® and (b) N = 10®. Isobars for 
larger N are similar to those in (b). 


E. Isothermal Compressibility, Thermal Expansion 
Coefficient and Heat Capacity at Constant Pressure 


In dimensionless reduced units Eq. (54) becomes 




p{x,N) 


dp{x,N) 

dx 


(56a) 


where the reduced isothermal compressibility kx is 


kt = 



(56b) 


One also has 



xp. 


ktp = ktp = 


(57) 
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The ideal gas exhibits 


ktP = 1 (ideal gas), (58) 


to which Kj'p for the Bose gas must asymptote for x —>• 
oo. 

We now derive an expression for dp (f) in terms of quan¬ 
tities already calculated. Writing p = at con¬ 

stant pressure one has the differential 


dp = 0 = (59) 


dt 


dj 


yielding 



p 


dp{t,'y,N) 

dt 

dp{t,'y,N) 

O7 


(60) 


Then using Eqs. (50c), (55b) and (60) one obtains 


^ d[xp{x,N)] 

7 \7^/ 


(61) 


where kt/ 7^ is given in Eq. (56a). Also, one has 


apT = X — . (62) 

7 

The ideal gas shows UpT = 1, to which x{ap/j) for the 
Bose gas must approach for a; —>■ 00. 

Finally, the difference Cp — Cy between the heat capac¬ 
ities at constant pressure and constant volume satisfies^^ 

Cp-Cy = -(63) 

/^x 

In dimensionless reduced parameters one obtains 




Cp - Cy (dp/7)^ 

NkB (« t / 7 ")' 


(64) 


For the ideal gas this quantity equals unity, which the 
Bose gas must approach for x ^ 00 . 

Shown in Fig. 12 are plots of kt/ 7^, dp/7 and (Cp — 
Cy)/NkB versus a; for (a) N = 10^ and (b) N = 10^. All 
three quantities show unphysical divergences and then 
negative values as x decreases into the BEC regime (not 
shown). Also shown are the associated plots of ktp and 
apT from Eqs. (57) and (62), respectively. One sees that 
(Cp — Cy)/[NkB), ktP and apT approach the same re¬ 
spective ideal gas value of unity at large a:, as required. 


V. RESULTS: CANONICAL ENSEMBLE 

In previous sections we pointed out a number of un¬ 
physical or unexpected predictions of the GCE formal¬ 
ism when N and x are both small (in the BEC regime) 
in addition to the known unphysically large fluctuations 
in N at small x even in the thermodynamic limit. In 


FIG. 12: (Color online) Plots of reduced isothermal compress¬ 
ibility thermal expansion coefficient Op/y, and the dif¬ 

ference {Cp — Cy)/NkB between the heat capacity at constant 
pressure and at constant volume for y = 1 and boson number 
(a) N = 10® and (b) N = 10®. Also shown in each panel are 
the products a^T and ktP which for an ideal gas are both 
equal to unity, as verified in the respective high-* limits in 
(a) and (b). The figure legend in (a) also applies to (b). 


this section we resolve these problems by calculating the 
thermodynamics using the CE formalism which can give 
exact results for a finite system with fixed N in thermal 
contact with a temperature reservoir. The partition func¬ 
tion Q{N) and average number of bosons fii in a given 
quantum state with energy Ei are calculated recursively 
as described in Sec. IIC, and our calculations are carried 
out with a maximum boson number N = 1000. Some of 
the thermodynamic properties for N = 1000 will be com¬ 
pared with the above unphysical and/or incorrect results 
predicted by the GCE formalism. 
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FIG. 13: (Color online) Fractional occupations Nq/N of the 
ground state versus x = -yt for four N values as determined 
from the CE formalism using Eq. (31b) (solid curves). The 
data for N = 1000 from Fig. 3 calculated using the GCE 
formalism are shown for comparison (dashed curve). 


A. Population Statistics 

The fractional occupancies Nq/N of the ground state 
for = 1, 10, 100 and 1000 obtained using Eq. (31b) 
are plotted versus x = yt in Fig. 13 (solid curves). On 
comparing the data with those in Fig. 3 for larger N 
values, one sees that the crossover between weak and 
strong increases of Nq /N versus x at x = 1 becomes much 
less well defined for small N. The data for N = 1000 
from Fig. 3 obtained from the GCE formalism are shown 
as the dashed curve in Fig. 13 for comparison. One sees 
that the CE and GCE formalisms are in reasonably good 
agreement for this value of N. 

The ratios Ni/N^ of the populations of a quantum 
state in each of the first four excited energy levels i = 1—4 
to that in the ground state Nq are plotted versus x in 
Fig. 14. Compared with the larger-A data in Fig. 4, 
the excited state populations for small N approach the 
ground state population at much larger x values than for 
larger N. 


B. Helmholtz Free Energy, Entropy and Internal 
Energy 

Within the CE formalism, we use the same definitions 
of t and X as given above in Eqs. (6) and (8) and of 
the ratio E/kNT in Eq. (9), respectively. To simplify 
notation we also define 

InQ(x) = (65) 

where Q{x) is calculated as described previously in 
Sec. lie. The reduced Helmholtz free energy F is given 


by 


and the reduced entropy S by 


S{x) 


S 

Nkn 


dF 


In Q{x) + X 


din Q(x) 
dx 


(67) 


The exact entropy at t —?► 0 with fixed y is easily ob¬ 
tained for any finite N. For x —>■ 0, only the ground state 
term with Ux = ny = 1 in Eq. (30c) is signiheant. Fur¬ 
thermore, when calculating Q{N) we must hold both T 
and A constant for each N. Then the factor Qi{k) is 


Qi(k, X —>■ 0) = exp 



(68a) 


where the expression for g is given in Eqs. (10). Inserting 
this into Eq. (30a) and carrying out the sum over k yields 


Q{N,: 


lnQ{N, X —>■ 0) = 


0) = exp 
hiQ 


r 2A1 


r 2a(TV)1 

. 9 . 

= exp 

X 


2a(fV) 


(68b) 


N Nx ■ 

Then the reduced free energy is obtained from Eq. (66) 
as 


F(x = 0) = 


2a{N) 
Ny ’ 


(69) 


where we used the dehnition x = yt. From the relation 
S = —dF{t,y,N)/dt one obtains the zero-temperature 
entropy 


S{t = 0) = 0, (70) 


which is valid for arbitrary finite N. This result makes 
physical sense, because there is only one way to put N 
indistinguishable bosons into an orbitally nondegenerate 
ground state. 

The reduced entropy within the CE formalism ob¬ 
tained from Eq. (67) is plotted versus x at small x < 1 
for TV = 1, 10, 100 and 1000 in Fig. 15. Also shown as 
the dashed curve are the data for TV = 1000 obtained 
from the GCE formalism in Fig. 6(a). One sees that the 
(incorrect) finite value of S' for T —>■ 0 and TV = 1000 
obtained with the GCE formalism is corrected using the 
CE formalism. 

The reduced internal energy U is 


or 


U = 


U 

NknT-E 


F + tS = t^y 


dlnQ 
dx ’ 


(71) 


U dlnQ 

Nk^T dx 


(72) 


where we used the relation x = yt. We find C'v(T —>■ 
0) = 0 and do not present plots of Cy{t) because they 
are similar to those in Fig. 5 obtained using the GCE 
formalism. 
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FIG. 14: (Color online) The ratio Nq/N versus x from Fig. 13 and the ratios Ni/No of the number of bosons Ni occupying a 
state in each of the first four excited energy levels with quantum numbers n^+riy = 5, 8, 10 and 13 in Eq. (7) to the ground-state 
occupation number No for (a) = 1, (b) N = 10, (c) N — 100 and (d) N = 1000. For the ground state + riy = 2. 


C. Pressure 

The reduced pressure p within the CE formalism is 
given by 


- _ _P}^_ _ _ f 

^ IV/cbTe V ^7 /1 dx 

The compression factor is 


(73) 


- pA _7 dhiQ{x) 
^ ^ NkBT ^ di 


(74) 


Comparing Eqs. (74) and (72) demonstrates that 



(75) 


which says that the pressure is equal to the average en¬ 
ergy density. This type of relationship is expected from 


dimensional considerations. For 3D Bose and Fermi gases 
in the thermodynamic limit, one obtains the similar ex¬ 
pression p = {2/3)U/V, where V is the volume of the 
gas."^ Plots oi p = versus x within the CE formal¬ 

ism are similar to those of the GCE formalism in Fig. 7(a) 
and 7(b) and are therefore not presented here. 

We solve Eq. (73) parametrically using x as an implicit 
parameter. We first calculate InQ(x). Then at constant 
7, one has t = xjj for a pressure versus temperature 
isochore, whereas at constant t one has 7 = x/t for a 
pressure versus area isotherm. Shown in Fig. 16 are p 
versus 7 isotherms at t = 0.5, 1 and 1.5 for = 1, 10, 
100 and 1000 in panels (a), (b), (c) and (d), respectively. 
As N increases, a hump appears at the crossover area 
7 ~ 1 for = 100 that is clearly defined by A^ = 1000. 
An important feature of these plots is that the slope is 
always negative. This means that kt is always finite 
and positive. This behavior is in contrast to the data 
for N = 1000 in Fig. 16(d) obtained using the GCE for- 
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FIG. 15: (Color online) Normalized entropy per boson S/NkB 
versus x = yt in the small-* regime for = 1 to 1000 ob¬ 
tained within the canonical ensemble (CE) formalism using 
Eq. (67) (solid curves). Also shown is the incorrect predic¬ 
tion for N = 1000 from Fig. 9(a) obtained within the grand 
canonical ensemble (GCE) formalism (dashed curve). 


malism, where one sees maxima in p versus 7 at 7 ~ 1 , 
which causes kt to exhibit an unphysical divergence on 
reducing 7 towards 7 ~ 1 and then unphysical negative 
values at lower 7 values. 


TABLE III: Reduced pressure p for reduced areas 7 = 0.5, 
1 and 1.5 and reduced isothermal compressibility all 

at zero temperature, versus N as predicted by the cononical 
ensemble formalism via Eqs. (78) and (82), respectively. 


login N 

p{t ^ 0) 

7 - 0.5 

p{t = 0) 

7^1 

p(t = 0 ) 

7 — 1.5 

KTfi = 0)/7^ 

0 

3.3231E-I-00 

8.3077E-01 

3.6923E-01 

6.0185E-01 

1 

8.8900E-01 

2.2225E-01 

9.8778E-02 

2.2497E+00 

2 

1.7417E-01 

4.3542E-02 

1.9352E-02 

1.1483E+01 

3 

2.S031E-02 

7.0078E-03 

3.1146E-03 

7.1349E-I-01 

4 

3.9907E-03 

9.9768E-04 

4.4341E-04 

5.0116E+02 

5 

5.2505E-04 

1.3126E-04 

5.8339E-05 

3.8092E-I-03 

6 

6.5526E-05 

1.6382E-05 

7.2807E-06 

3.0522E+04 

7 

7.SS14E-06 

1.9704E-06 

8.7571E-07 

2.5376E-I-05 

8 

9.2284E-07 

2.3071E-07 

1.0254E-07 

2.1672E+06 

9 

1.0589E-07 

2.6471E-08 

1.1765E-08 

1.8888E-I-07 

10 

1.1959E-08 

2.9897E-09 

1.3287E-09 

1.6724E+08 

11 

1.3353E-09 

3.3383E-10 

1.4837E-10 

1.4978E-I-09 

12 

1.4722E-10 

3.6804E-11 

1.6357E-11 

1.3585E+10 

13 

1.6115E-11 

4.0289E-12 

1.7906E-12 

1.2410E+11 

14 

1.7499E-12 

4.3747E-13 

1.9443E-13 

1.1429E+12 

15 

1.8905E-13 

4.7262E-14 

2.1006E-14 

1.0579E+13 

16 

2.0311E-14 

5.0777E-15 

2.2568E-15 

9.8470E+13 

17 

2.1666E-15 

5.4165E-16 

2.4073E-16 

9.2311E+14 


As noted in the introduction, a nonzero pressure must 
occur at t = 0 in a noninteracting Bose gas in a 2D box 
with Dirichlet boundary conditions because the ground 
state energy depends on the area.^^ At t = 0, all N 


bosons are in the ground state with = riy = 1. From 
Eq. (5a), the energy of the ground state {rix = Uy = 1) 
containing N bosons at t = 0 is 


NEo{t = 0) _ 2a{N) _ 2a{N)AE 
^bTe 7 ^ 

Then the pressure p at t = 0 is 

p{t = 0) d{NEo/kBTE) 2a{N)AE 

/cbTe dA A2 


(76) 


(77) 


Using the definitions 7 = {A/N)/{A/N)e = AjAE and 
of p in Eq. (73), one obtains the reduced pressure 


p{t = 0 ) 


2a{N) 
Ny^ ■ 


(78) 


Values of p{t = 0) for N = 10° to 10^^ obtained from 
Eq. (78) using the values of a{N) in Table I are listed in 
Table III for 7 = 0.5, 1 and 1.5. One sees that p{t = 0) 
is quite large for small N, but decreases rapidly as N 
increases. Isochores of p versus t for 7 = 0.5, 1 and 1.5 
obtained using Eq. (73) are shown at low t for N = 1, 10, 
100 and 1000 in panels (a), (b), (c) and (d) of Fig. 17, 
respectively. One indeed sees that p{t -y 0) decreases 
rapidly with increasing iV, with the t = 0 values in 
agreement with those listed in Table III. Also shown 
in Fig. 17(d) are the corresponding isochores in Fig. 9(a) 
obtained from the GCE formalism (dashed curves), for 
which the incorrect limit p{t —?► 0 ) = 0 is obtained for 
each of the three 7 values. 


D. Isothermal Compressibility, Thermal Expansion 
Coefficient, Heat Capacity at Constant Pressure 


The reduced isothermal compressibility kt /y'^ is 


'i-i , 

-^[x) = kt 
7 


Nk^TE 

y^VE 



1 


3 d^InQ(x) ’ 


(79) 

where the large-* ideal gas limit ktP = ktP^x -y 00 ) = 1 
is expected. 

The reduced thermal expansion coefficient cxp/y is 


oip, . OpTs 1 

— (*) = —— =- y , 

7 7 7^ 


o d\nQ 
dx 

dx'-^ 


1 

5 

X 


(80) 


where the large-* limit is expected to be the ideal gas 
value OfpT = dp*)* -y 00 ) = 1. The normalized differ¬ 
ence between the heat capacities at constant pressure Cp 
and constant volume (constant area) Cy is given by 


Cp-Cy 

NkE 


(*) 


kt/7^ 


(81) 


All of these quantities are plotted versus * in Fig. 18 for 
N = 1, 10, 100 and 1000 in panels (a)-(d), respectively. 
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FIG. 16: (Color online) Reduced pressure p = pyE/ksTs versus reduced area 7 = v/ve isotherms at reduced temperatures 
t = T/Te = 0.5, 1 and 1.5 and boson numbers (a) = 1, (b) N = 10, (c) N = 100 and (d) N = 1000 (solid curves) obtained 

using the canonical ensemble (CE) formalism. Also shown in (d) are the corresponding incorrect data obtained using the grand 
canonical ensemble (GCE) formalism for N = 1000 in Fig. 10(a) (dashed curves). Note that the scale of the ordinate in (d) is 
different than in (a)-(c). 


One sees that with the CE formalism, one does not en¬ 
counter the unphysical divergences and other inaccura¬ 
cies discussed above that occur with the GCE formalism 
at small x and N values where BEC comes significant. 

Using Eq. ( 68 b) for \iiQ{x —>■ 0) together with the 
general definition for kt/ 7 ^ in Eq. (79), one obtains the 
zero-temperature limit at fixed 7 given by 


7 


N 

4a(iV) ■ 


(82) 


A list of values of 0) versus N obtained using 

Eq. (82) is given in the fifth column of Table III, where 
to calculate these we used the a{N) values in Table I. 
Similarly, we find that (Q;p/ 7 )(t = 0) = 0 and hence (Cp — 
Cy){t = 0) = 0 using Eq. (81). These zero-temperature 


results are in agreement with the a: = 0 limits of the 
respective plots in Fig. 18. 


VI. SUMMARY 

We confirmed the literature result that BEC does not 
occur in the thermodynamic limit at finite temperature T 
in a noninteracting Bose gas confined to a 2D box. How¬ 
ever, as also previously reported for finite N, BEC does 
occur in 2D, where the ground state boson occupation 
is Aq/YV —1 at fixed area A for T —>■ 0, but without 
any phase transition occurring.^® The lack of a phase 
transition is confirmed from the analytic behavior of the 
calculated Cv(T) upon traversing the characteristic tern- 


















19 


LU 



LU 

> 

a. 




t 


t 




t 


t 


FIG. 17: (Color online) Reduced pressure p = pvY^/{k^TY) versus reduced temperature t isochores at low t at fixed reduced 
areas 7 = 0.5, 1 and 1.5 and boson numbers (a) = 1, (b) N = 10, (c) N = 100 and (d) N = 1000 (solid curves) obtained 

using the canonical ensemble (CE) formalism. Also shown in (d) are the corresponding data obtained using the grand canonical 
ensemble (GCE) formalism (dashed curves) from Fig. 9(a). The axis scales are different in each plot in order to emphasize the 
nonzero values of p{t 0) obtained with the CE formalism. 


perature Te. Thus the parameter Te that we define cor¬ 
responds to a crossover temperature between weak and 
strong increases in Nq/N and in the low-lying excited 
states with decreasing T at fixed A and not to a phase 
transition temperature. We find that Te decreases with 
increasing N according to Te ^ l/log(fV) at fixed area 
per boson (A/7V)e yielding Te{N —>• 00 ) = 0. Hence 
BEG is precluded at finite T in the thermodynamic limit 
in 2D whereas it does occur at low T with finite N in the 
absence of a BEG phase transition, a perhaps counterin¬ 
tuitive result. 

The main contribution of this paper is a comprehensive 
and detailed study of the thermodynamic properties of 
noninteracting bosons in a 2D box with Dirichlet bound¬ 
ary conditions. Such a study has not been carried out 
before to our knowledge and is therefore a benchmark 


for future studies on similar systems. We used both the 
GGE and GE formalisms for the calculations. The GGE 
formalism generally gives accurate results for the ther¬ 
modynamic properties at large N and large values of the 
product TA, but fails to give correct results for small N 
at small TA values where significant BEG occurs. Such 
failures of the GGE formalism in the latter ranges of pa¬ 
rameters include incorrect predictions of nonzero entropy 
and zero pressure, strong deviations of the ratio of the 
pressure to the energy density pj(UjA) from the exact 
GE value of unity, and divergent and/or negative values 
of kt, ttp and Cp. These incorrect behaviors predicted by 
the GGE formalism are revealed using the GE formalism 
which permits numerically and analytically exact results 
to be obtained, albeit at comparatively small N. Thus 
apart from the specific study reported here, we hope that 
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FIG. 18: (Color online) Reduced isothermal compressibility kt/ 7 ^, thermal expansion coefficient dp/ 7 , and normalized differ¬ 
ence (Cp — Cy)/NkB between the heat capacity at constant pressure and at constant volume versus x = 'yt for boson numbers 
(a) N = 1, (b) N = 10, (c) N = 100 and (d) N = 1000. Also shown in each panel are the products apT and ktP which for 
an ideal gas are both equal to unity, as expected and found in the respective large-x limits in (a)“(d). The figure legend in (a) 
applies to all panels. Note the different ordinate scales for each panel. 


the present results will be more generally useful because 
they illustrate several generic shortcomings of the GCE 
formalism in predicting the thermodynamic properties of 
finite quantum boson systems. 
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